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Dynamics and structure of 
turbulent premixed flames 

By R. W. Bilger 1 , N. Swaminathan 1 , G. R. Ruetsch, AND N. S. A. Smith 


1. Motivation and objectives 

In earlier work (Mantel &; Bilger, 1994) the structure of the turbulent premixed 
flame was investigated using statistics based on conditional averaging with the re- 
action progress variable as the conditioning variable. The DNS data base of Trouve 
and Poinsot (1994) was used in this investigation. Attention was focused on the 
conditional dissipation and conditional axial velocity in the flame with a view to 
modeling these quantities for use in the conditional moment closure (CMC) ap- 
proach to analysis of kinetics in premixed flames (Bilger, 1993). Two remarkable 
findings were made: there was almost no acceleration of the axial velocity in the 
flame front itself; and the conditional scalar dissipation remained as high, or higher, 
than that found in laminar premixed flames. The first finding was surprising since 
in laminar flames all the fluid acceleration occurs through the flame front, and this 
could be expected also for turbulent premixed flames at the flamelet limit. The find- 
ing gave hope of inventing a new approach to the dynamics of turbulent premixed 
flames through use of rapid distortion theory or an unsteady Bernoulli equation. 
This could lead to a new second order closure for turbulent premixed flames. The 
second finding was contrary to our measurements with laser diagnostics in lean 
hydrocarbon flames where it is found that conditional scalar dissipation drops dra- 
matically below that for laminar flamelets when the turbulence intensity becomes 
high. Such behavior was not explainable with a one-step kinetic model, even at 
non- unity Lewis number. It could be due to depletion of H2 from the reaction zone 
by preferential diffusion. The capacity of the flame to generate radicals is critically 
dependent on the levels of H2 present (Bilger, et al . , 1991). It seemed that a DNS 
computation with a multistep reduced mechanism would be worthwhile if a way 
could be found to make this feasible. 

Truly innovative approaches to complex problems often come only when there is 
the opportunity to work close at hand with the (in this case numerical) experimental 
data. Not only can one spot patterns and relationships in the data which could be 
important, but one can also get to know the limitations of the technique being 
used, so that when the next experiment is being designed it will address resolvable 
questions. A three-year grant from the Australian Research Council has enabled 
us to develop a small capability at the University of Sydney to work on DNS of 
turbulent reacting flow, and to analyze data bases generated at CTR. Collaboration 
between the University of Sydney and CTR is essential to this project and finding 
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a workable modus operandum for this collaboration, given the constraints involved, 
has been a major objective of the past year’s effort. 

The overall objectives of the project are: (1) to obtain a quantitative understand- 
ing of the dynamics of turbulent premixed flames at high turbulence levels with a 
view to developing improved second order closure models; and (2) to carry out new 
DNS experiments on turbulent premixed flames using a carefully chosen multistep 
reduced mechanism for the chemical kinetics, with a view to elucidating the laser 
diagnostic findings that are contrary to the findings for DNS using one-step ki- 
netics. In this first year the objectives have been to make the existing CTR data 
base more accessible to coworkers at the University of Sydney, to make progress on 
understanding the dynamics of the flame in this existing CTR data base, and to 
carefully construct a suitable multistep reduced mechanism for use in a new set of 
DNS experiments on turbulent premixed flames. 

2. Accomplishments 

2.1 Accessing the data base 

A Fortran 77 program has been written that allows easy access to the DNS data 
base of Trouve and Poinsot (1994) from the computers at NAS and also allows 
efficient computation of derivatives of the data and conditional statistics. It can 
also be used on the DEC Alpha workstations used at the University of Sydney 
for this work, but with the limitation that data transfer and storage is limited 
and further coding has to be written for converting the subroutines that calculate 
derivatives. So far we have only one field at one time for the unity Lewis number 
case available in Sydney at 32-bit accuracy. Some processing has been achieved 
with this at Sydney, but more memory is needed before we can do processing that 
requires several arrays to be stored at the one time. A further 64 MBytes of RAM 
will become available soon and this will remove this limitation. It is proposed to 
write most of the Trouve and Poinsot data base on tape at 32 bit accuracy and 
transport it to Sydney in this form. 

2.2 Dynamics of the flame 


2.2.1 Introduction 

During the time that the senior research fellow (Bilger) was at CTR, a good deal 
of progress was made investigating the velocity and pressure field in the Trouve and 
Poinsot flame. It was found that the pressure difference across this flame is small 
and the acceleration of the mean flow comes from the normal Reynolds stress. This 
dominance of the normal Reynolds stress makes it seem likely that the use of rapid 
distortion theory or an unsteady Bernoulli equation will not be successful. 

Efforts to progress further were frustrated by apparent anomalies in any balances 
that contained the reaction rate. This problem took some time to uncover as most 
balances involve time dependent terms, and these were unavailable without running 
the original code, the new Fortran 90 version not then being available. The flame is 
not statistically stationary and the time dependent term in any averaged equation 
can be quite significant. An equation for the dilatation in the flow was derived which 
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has no unsteady term, and this showed that there was an error in the parameters 
being used. This has only been cleared up in the last month or so. The limitations 
on accessing the data base at Sydney are such that limited processing has only been 
possible up to now, but it is expected that this problem will be overcome soon. 

The results for the dilatation equation are interesting in themselves and are pre- 
sented here in detail. 

2.2.2 Dilatation equation 

Instantaneous reaction rate in turbulent flows is related to local value of scalar 
dissipation rate as demonstrated by Bilger (1976) in the fast chemistry limit and 
by Peters (1983) for finite rate chemistry for diffusion flames. In the high Da 
(Damkohler number) limit, Bray (1980) showed the same dependency for premixed 
flames. Hence, scalar dissipation rate is a crucial quantity in reacting flows. Re- 
cently, this quantity has attracted a quite a few modeling studies (Chen et al. 1989, 
Girimaji 1992). Mantel and Bilger (1994) analyzed DNS data base of Trouve and 
Poinsot (1994) to comprehend the behavior of scalar dissipation rate in turbulent 
premixed flames. They observed that the behavior of scalar dissipation rate is in- 
dependent of the position inside the turbulent flame brush. Here, we demonstrate 
that the conditional scalar dissipation rate can be obtained from the conditional 
dilatation equation for premixed flames, knowing the probability density function 
of the progress variable. 

To ease the comparisons of the results with DNS, we make all the quantities in 
the following discussion dimensionless, unless otherwise specified, by using acous- 
tic scales as in Trouve and Poinsot (1994). One can write the mass conservation 
equation as 

V-Un^ + pU-Vfl/p), (1) 

where p and u respectively denote density and velocity vector. By making use of 
the equation of state and defining the reaction progress variable c as 

c = (T-T u )/(T b -T u ) = — ®[( 7 -l)T-l], 

ol 

where a is a heat release parameter (Williams 1985), one can write Eq. (1) as 


V 


u = 


a 

1 — a 



+ pu • Vc . 


( 2 ) 


Substituting the governing equation for c in Eq. (2), one can obtain 


Vu — 


a 


1 -c* 


“ e+ ^ v '‘ vc 


(3) 


where Cj c and p are respectively the reaction rate of c and the local dynamic viscosity 
of the fluid (Trouve & Poinsot 1994). This equation is referred to as the dilatation 
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Figure 1. Variation of average velocity (il), reaction rate (u; c ), and diffusive flux 

along the inhomogeneous direction x. u: , DNS; , Eq.(4), unmarked 

line diffusive flux. See equation (4). 


equation in the following discussions. Averaging this equation, after noting the 
statistical one-dimensionality of the DNS (Trouve &; Poinsot 1994), yields 

u{x,t) = u(-oo,t) + J^u c (x' ,t)dx' + ^ 

The average velocity u(x,t) obtained from the above equation is compared with 
DNS results in Fig. 1 for Re = 1000, Pr = 0.75, Le = 1, and a = 0.75 at t = 4.5. 
Time is dimensionless with respect to initial eddy turnover time (Trouve &: Poinsot 
1994). The agreement is excellent and encourages us to proceed further. 

2.2.3 Conditional dilatation equation 

Following Klimenko (1990), the conditioning process can be expressed using a 
Dirac delta function # = 6(c[x,t\ - £), where c and £ are respectively progress 
variable and its sample space variable. Using this notation, the conditional average 
of any quantity B can be written as < > = < B\( > where denote 

the probability density function of progress variable. After applying the rules for 
the differentiation of the delta function, we can obtain the conditional dilatation 
equation as 


V . [< u|C > P<] = - (< u • Vc lC > Pc) + > p < 


+ 


a 


1 


( 5 ) 


1 — a RePr 


D<P<, 
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where D $ is the conditional diffusion expressed as =< V • /xVc|£ >. For statis- 
tically 1-D flow (as in the DNS) Eq. (5) can be written as 


= RePr 


1 — a 


a 


^(<t 1 |C>P<)+^(<uVc|C>P { ) 


- RePr < u; c |C > P( 


( 6 ) 


The second part of the above equation is obtained by relating the conditional dif- 
fusion to conditional dissipation N^. In principle, one can get the conditional dis- 
sipation rate from the above equation. 

2.8 Multistep reduced mechanism 


2.3.1 Formulation 

Possibilities of direct numerical predictions of turbulent reacting flows, even in 
simple geometries, with multistep elementary kinetic mechanisms are remote with 
the current computational hardware. Hence, researchers often simplify the chemical 
reactions to a single global step. Direct simulations with single-step chemistry 
(Trouve Poinsot 1994, Swaminathan et a/., 1995, Givi 1989), although simplified 
in certain sense, have given us valuable insight into the different physical processes 
involved in reacting flows. To further our understanding, direct simulations with a 
systematically reduced kinetic scheme would be of great interest. Mahalingam et 
at. (1995) have simulated turbulent nonpremixed flames using a two-step chemical 
scheme which is similar to the Zeldovich-Linan mechanism. Here, we present a 
systematically reduced two-step scheme (Peters & Williams 1987, Williams 1991) 
for hydrocarbon flames from the direct numerical simulation point of view. 

Peters (1985) has shown the strategy of reducing full kinetics to a simplified four- 
step mechanism. Simplification strategy consists of order of magnitude arguments, 
steady state, and partial equilibrium approximations for appropriate minor species. 
This reduced mechanism has been used in computational (Bilger et al. 1990, Pe- 
ters & Kee 1987) and asymptotic (Peters & Williams 1987, Seshadri & Peters 1988) 
studies of laminar flames. These studies improved our understanding of flame struc- 
ture and extinction mechanisms. For example, a laminar diffusion flame modeled 
by a single global step extinguishes at high strain rates after allowing fuel to leak 
through the reaction zone, whereas studies with four-step chemistry indicate that 
the oxidizer leaks through the reaction zone due to radical depletion on its rich side. 
These differences indicate that direct simulations with a reduced mechanism would 
further improve our understanding of turbulent flames and thereby allow us to con- 
struct more accurate models for engineering predictions. Due to the computational 
requirements, we further reduce a four-step mechanism to two steps as discussed 
below. Reducing a full mechanism to four steps can be found elsewhere (Peters 
1991, Bilger et al. 1990). 
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Table 1. Elementary steps involved and their rate constants 


Step 

Reaction 

A 

b 

E 

1 

CH\ 4" H — ♦ CH 3 + H 2 

6.6E+08 

1.6 

10840.0 

2 

CO + OH ^ C0 2 + H 

1.2E+07 

1.4 

-730.0 

3 

H + 0 2 + M —+ H0 2 + M 

3.6E+17 

-0.72 

0.0 

4 

H + 0 2 ^ OH + 0 

8.3E+13 

0.0 

14413.0 


K = AT b exp[— E/RT], cal-mol-cm-sec-K 


Here we start with the four steps given by Bilger et al (1990). These four steps 
are 


CHi + 2H + H 2 0 — » CO + 4 H 2 

I 

CO + H 2 0 ^ C0 2 4- H 2 

II 

2 H 2 + 0 2 —+ 2 H 2 0 

III 

ZH 2 + 0 2 — 2 H 2 0 + 2 H. 

IV 


The rates of these four steps are given as linear combinations of some elementary 
reactions. Excluding all but the most elementary steps, the rate expression given 
by Bilger et al (1990) can be expressed as 


W/ = i, 


VJJ = U>2, 


VI J I = u> 3 , 


VIV — V 4 , 


where the roman and arabic subscripts respectively denote global and elementary 
steps. These elementary steps with their rate constants are given in Table 1. 

The pressure dependence of elementary step 3 allows us to make a steady state 
approximation for H atom at pressures typically above one atmosphere. This as- 
sumption renders 

2 K [MM 

c {H 2 0} 2 


1 - A 


[CH 4 

[ 02 ] 


(7) 


where K c is a combination of equilibrium constants of elementary reactions involved 
in making partial equilibrium approximation of OH and steady state approximation 
for O atoms (Bilger et al. 1990, Peters & Kee 1987). The ratio of rate constants 
of elementary steps one to four is denoted by A and its magnitude is about four 
for a temperature range of interest. Stoichiometry of the resulting reaction rate 
expressions of individual species lead us to a three-step mechanism given by (Peters 
& Williams 1987) 

Cff 4 + 0 2 — ♦ CO + H 2 + h 2 o 
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CO + h 2 o ^ co 2 + h 2 

2H 2 + 0 2 — » 2 H 2 0. 

This mechanism can further be reduced to two steps by making partial equilibrium 
approximations for the water-gas shift reaction. With this assumption, the reaction 
rates of individual species are given as 

U>CH 4 = — W/, WOj = — (w/ +w///) 

«« 2 = -2 (w/ - w///)/( 1 + a), u>co = -2a(d>/ - a>//j)/(l + a) 

u>co 2 = -[(1 - <*)<^/ + 2aw///]/(l + a), uh 7 o = — 2[aw/ + w///]/( 1 + a), 

where a is ratio of CO to H 2 concentrations at partial equilibrium of water-gas shift 
reaction and is assumed to be a constant across the reaction zone. Stoichiometry 
coefficients in the above rate expressions give the two steps as 

CH a + 0 2 — > 2 + ocCO) + -+C0 2 + ff 2 0 

1 + a 1 + a 1 + a 

O 2 + tt~[H 2 + aCO] — ► + - ^—h 2 o. 

1 + a 1 + a 1+a 

Laminar flame calculations (Bilger et al 1991) with a skeletal mechanism suggest 
that a is 0(1). Hence the above two steps become 

CH a + 0 2 — ► [J5T 2 + CO] + H 2 0 

0 2 + [H 2 + CO] — ► ff 2 0 + C0 2 . 

In this mechanism, [ff 2 + 00] clearly plays a role of intermediates. Hence, by 
denoting them as /, CH 4 as F, 0 2 as Oxi , H as iZ, and the remaining species as 
product P, we get 

F + Oxi — ► I + P f 

Oxi + I —> 2 P. //' 

Reaction rates of these steps, according to the elementary reactions involved, are 
given as uy = K\ [F][P] and u n * = K' 3 [0xi][R}. Radical concentration [R] is given 
by a modified form of Eq. (7) as suggested by Peters (1995): 

[R] = AV[0xz]°* 5 [/] 1 * 5 exp(-aA[P]/[0xi]), 

where K r is related to the equilibrium constant K c , and a is a constant. This form 
is used to avoid the discontinuity which is implicit in the steady state approximation 
for H atom (see Eq. (7)) as this will give problems in DNS. By matching the flame 
speed eigen value for premixed flames (Peters & Williams 1987) or extinction scalar 
dissipation rate at stoichiometric mixture fraction for diffusion flames (Seshadri & 
Peters, 1988), one can obtain a = -^15/4. One can also show that the first and 
second step of the above two-step mechanism release 40% and 60%, respectively, of 
the overall heat release. Structure of nonpremixed flames with the above two steps 
is presented in the following discussions. 
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FIGURE 2. Structure of a laminar (Tsuji) methane-air diffusion flame for a strain 
rate value of 100 sec" 1 . 



Z 

FIGURE 3. Structure of a laminar diffusion flame used as an initial field for DNS 
of turbulent nonpremixed flame. 

2.S.2 Nonpremixed flame structure 

Calculation of Tsuji type laminar diffusion flames axe carried out to understand 
the laminar flame structure and its relation to different rate constants involved. 
Figure 2 depicts the structure of a methane-air flame at a strain rate of 100 sec~ l 
in mixture fraction space. 

The rate constants are derived from those given in Table 1. The specific moles are 
normalized by the free stream value of fuel species. Radical concentration is about 
ten times lower than the intermediate concentration. The intermediate is formed 
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FIGURE 4. Turbulent nonpremixed flame structure with two-step kinetics at 
t = 1.5. Time t is made dimensionless with eddy turnover time at scalar field 
initialization. 


on the rich side of the reaction zone while it is consumed on the lean side. The 
maximum temperature is about 1900 K, while the free stream value is 300 K. The 
size of fuel depletion zone in mixture fraction space is about one fifth of oxidizer 
consumption zone. This flame structure is consistent with Bilger et aV s (1991) 
calculation using a skeletal mechanism. 

The value of stoichiometric mixture fraction can be increased, to ease resolution 
requirements in DNS, by diluting the reactant streams. The combination of an 
oxidizer stream having 50% 0 2 and 50% iV 2 by weight and a fuel stream of 25% 
methane and 75% Argon by weight has a mixture fraction stoichiometric value of 
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1/3. These mixtures have equal moleculax weights and densities which are the same 
as that of air. These attributes are attractive for direct simulations of turbulent 
nonpremixed flames with low heat release. The structure of a laminar diffusion flame 
for the above reactants is shown in Fig. 3. The rate constants for this calculation 
are evaluated at 1800 K. This structure is used to initialize DNS calculations using a 
psuedospectral algorithm (Swaminathan et al. 1995). The structure of the turbulent 
flame in mixture fraction space is shown in Fig. 4 in the form of scatter plot. The 
scatter in the values of specific moles is due to unsteady effects. More analysis of 
this calculation is expected to guide us to design a better experiment. 

3. Future work 

Our immediate goal is to settle the parameter values to be used in the two-step 
reduced mechanism for DNS. Laminar and two-dimensional turbulent premixed 
flames will be tested with the new Fortran 90 compressible code (see Ruetsch and 
Broadwell 1995, Smith 1995), and then ‘production’ runs will be made in three 
dimensions with various values chosen for the Lewis number of the intermediate. 
The simulation conditions chosen will be similar to those of Trouve and Poinsot but 
hopefully with an inflow boundary condition of non-decaying turbulence. Given 
sufficient computational resources, it is hoped that the ‘production’ runs can be 
completed in time for analysis in the 1996 CTR Summer Program. 

We will also pursue the goal of making the whole of the existing data base available 
in Sydney. As soon as the expanded memory for our DEC Alpha is available we 
will return to analysis of the conditionally averaged momentum, energy, Reynolds 
stress, and Reynolds flux equations. It is still considered likely that a new approach 
to second order closure will be found. 
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